clear; clc;
blh=Simulation.BlochSimulator.BlochSimulator('Rp', 2*pi*1e3, 'OmegaY', 0, 'OmegaX', @(t) 2*pi*4.7e3*cos(2*pi*70e3*t), 'OmegaZ', 2*pi*70e3);
tic; [mVect, tList] = blh.solveTransient(0:1e-6:1e-3, [0; 0; 0.5]);toc
blh.plot();


%%
clear; clc
rb85 = Miscellaneous.GetRb85;
rb87 = Miscellaneous.GetRb87;

rfm85 = Simulation.BlochSimulator.RFMagnetometer(rb85);
rfm87 = Simulation.BlochSimulator.RFMagnetometer(rb87);

rfm85.setPumpBeam(100);
rfm87.setPumpBeam(100);

bList = linspace(0, 20e-6, 301);
dmdList85b = arrayfun(@(b0) calcDemod(rfm85, b0, 100e-9, 47e3), bList);
dmdList87b = arrayfun(@(b0) calcDemod(rfm87, b0, 100e-9, 47e3), bList);

fList = linspace(1e3, 100e3, 300);
dmdList85f = arrayfun(@(f) calcDemod(rfm85, 10e-6, 100e-9, f), fList);
dmdList87f = arrayfun(@(f) calcDemod(rfm87, 10e-6, 100e-9, f), fList);

figure;
subplot(2, 1, 1);
plot(bList, [dmdList85b.R]*rb85.parameter.abundance, 'ro-', ...
     bList, [dmdList87b.R]*rb87.parameter.abundance, 'bd-', ...
     bList, [dmdList85b.R]*rb85.parameter.abundance + [dmdList87b.R]*rb87.parameter.abundance, 'ks-')
grid on; xlim(minmax(bList));

subplot(2, 1, 2);
plot(fList, [dmdList85f.R]*rb85.parameter.abundance, 'ro-', ...
     fList, [dmdList87f.R]*rb87.parameter.abundance, 'bd-', ...
     fList, [dmdList85f.R]*rb85.parameter.abundance + [dmdList87f.R]*rb87.parameter.abundance, 'ks-')
grid on; xlim(minmax(fList));
%%
clear; clc
rb = Miscellaneous.GetRb85;
rfm = Simulation.BlochSimulator.RFMagnetometer(rb);
rfm.setPumpBeam(100).setMainField(10e-6).setDrivingField(500e-9, 47e3);
tic; rfm.solveTransient(0:1e-6:1000e-6, [0 0 0.5]);toc
rfm.plot()

%%
function demod = calcDemod(rfm, b0, b1, freq)
    rfm.setMainField( b0 )...
       .setDrivingField( b1, freq );
    rfm.solveSteadyState();
    demod = rfm.demodulation();
end


